
#pragma once

#include <thrust/fill.h>
#include <thrust/count.h>
#include <cusp/dia_matrix.h>

namespace cusp
{
namespace gallery
{

/** \brief Build the identity matrix.
*/
template <class MatrixType>
void eye(MatrixType& A, size_t num_rows, size_t num_cols,
	typename MatrixType::value_type value = typename MatrixType::value_type(1))
{

	typedef typename MatrixType::index_type IndexType;
	typedef typename MatrixType::value_type ValueType;

	//A.resize(num_rows, num_cols);

	size_t nnz = thrust::min(num_rows, num_cols);
	cusp::dia_matrix<IndexType, ValueType, typename MatrixType::memory_space> dia(num_rows, num_cols, nnz, 1);
	dia.diagonal_offsets[0] = 0;
	thrust::fill_n(dia.values.values.begin(), nnz, value);
	dia.num_entries = dia.values.values.size() - thrust::count(dia.values.values.begin(), dia.values.values.end(), ValueType(0));

	A = dia;
}


} // end namespace gallery
} // end namespace cusp

